Analysis date: 2020-01-10

Setup

Load libraries

library(plyr)
library(gtools)
library(openxlsx)
library(pheatmap)
library(reshape2)
library(progress)
library(Matrix)
library(Hmisc)
library(lemon)
library(ggpubr)
library(effsize)
library(ggbeeswarm)
library(ggfortify)
library(ggpmisc)
library(ggrepel)
library(readxl)
library(DESeq2)
library(TOSTER)
library(tidyverse)
library(vsn)
library(fdrtool)
library(limma)
library(apeglm)
library(IHW)
library(Rtsne)
library(biomartr)
library(biomaRt)
library(MultiAssayExperiment)
library(PMA)
library(gplots)
library(RColorBrewer)
library(grid)
library(ConsensusClusterPlus)
library(survival)
library(survminer)
library(cowplot)
library(matrixStats)
library(DEqMS)

Load data

source("/Volumes/sd17b003/Sophie/Analysis/Screen_analysis/Figure_layouts.R")
load("/Volumes/sd17b003/Sophie/Analysis/CLL_Proteomics/CLL_Proteomics_final/Proteomics_Git/Robjects/CLL_Proteomics_Setup.RData")
load("/Volumes/sd17b003/Sophie/Analysis/CLL_Proteomics/CLL_Proteomics_final/Proteomics_Git/Robjects/CLL_Proteomics_ConsensusClustering.RData")
colData(multiomics_MAE)$PG <- as.factor(CCP_group5[rownames(colData(multiomics_MAE))])
colData(multiomics_MAE)$CCP6_RNA <- as.factor(CCP_group6_RNA[rownames(colData(multiomics_MAE))])

Limma Proteomics and genetic alterations

Functions

Function with DEQMS

Calculate_Limma <- function(assay_data, mut){
  # prepare assay data
  assay_data[is.infinite((assay_data))] <- NA
  is.na(assay_data) <- NA

  stopifnot(mut %in% colnames(proteomics_tbl_meta_biomart))
  # prepare row data
  row_data <- left_join((rownames(assay_data) %>% enframe(value = "rowname")),
              (proteomics_tbl_meta_biomart %>% 
                dplyr::select(rowname, start_position, end_position, chromosome_name,  mean_position) %>% unique %>%
                 group_by(rowname) %>% 
                 dplyr::slice(1) %>% ungroup() ),
              by=c("rowname")) %>% dplyr::select(-name) %>% as.data.frame()
  rownames(row_data) <- row_data$rowname
  
  stopifnot(all(rownames(assay_data)==rownames(row_data) ))

  # prepare column data
  col_data <- left_join((colnames(assay_data) %>% enframe(value = "colname")),
            (proteomics_tbl_meta_biomart %>% dplyr::select(colname, mut) %>% unique() ),
            by=c("colname") ) %>% dplyr::select(-name) %>% as.data.frame()
  rownames(col_data) <- col_data$colname
  
  stopifnot(all(colnames(assay_data)==rownames(col_data) ))
  if( (length(unique(col_data[,mut])) ==2) & !is.logical(col_data[,mut]) & is.numeric(col_data[,mut]) ){
    col_data[,mut] <- as.logical(col_data[,mut] == max(col_data[,mut]) )
  }
  
  #Creating an expression set
  raw_dataE <- ExpressionSet(assayData = assay_data,
                             phenoData = AnnotatedDataFrame(col_data),
                             featureData = AnnotatedDataFrame(row_data))
  validObject(raw_dataE)
  if(!is.logical(pData(raw_dataE)[,colnames(pData(raw_dataE))==mut]))stop("The mutation you want to look at is not logical")
  raw_dataE
  
  # Limma
  limma_data <- raw_dataE
  comparison <- c("mut - wt")
  
  colnames(pData(limma_data))[colnames(pData(limma_data)) == mut] <- "condition"
  pData(limma_data) <- mutate(pData(limma_data),condition= if_else(condition, "mut","wt") )
  limma_data <- limma_data[, !is.na(pData(limma_data)$cond)]
  limma.cond <- factor(pData(limma_data)$condition, ordered = FALSE)
  
  contrast.matrix <- model.matrix( ~ 0 + limma.cond)
  colnames(contrast.matrix) <- gsub("limma.cond", "", colnames(contrast.matrix))
  
  limma.object <- eBayes(
    contrasts.fit(
      lmFit(limma_data, design = contrast.matrix),
      makeContrasts(contrasts = comparison, levels = contrast.matrix)
      )
    )  

  ##### DEQMS
  tmp <- metadata(multiomics_MAE)$protein_description %>% dplyr::select(`Gene Name`, MinPSMQuant) %>% 
    mutate(MinPSMQuant= as.numeric(MinPSMQuant))
  psm.count.table <- tmp[,2] %>% as.data.frame()
  rownames(psm.count.table) <- tmp$`Gene Name`
  
  limma.object$count = psm.count.table[rownames(limma.object$coefficients),"MinPSMQuant"]
  limma_DEQMS = spectraCounteBayes(limma.object)
  
  DEqMS.results = outputResult(limma_DEQMS,coef_col = 1)
    
  limma_results_i <- DEqMS.results
    limma_results_i <- subset(limma_results_i, !is.na(logFC))
    limma_results_i$comparison <- comparison
    limma_results_i$mut <- mut
    colnames(limma_results_i)[c(8,9,10, 14, 15, 16)] <- 
      c("limma.t_beforeDEQMS", "pvalue.limma_beforeDEQMS" , "fdr.limma_beforeDEQMS", "t", "pvalue.limma", "fdr.limma")
    limma_results_i$fdr <- limma_results_i$fdr.limma
    limma_results_i$pvalue <- limma_results_i$pvalue.limma
    return(limma_results_i)

}

Annotation function

Annotate_Limma_Results <- function(limma_results){
  
  fdr_hit_threshold <- 0.001
  fdr_candidate_threshold = 0.01
  fc_hit_threshold <- log2(1.5)
  fc_candidate_threshold <- log2(1.2)
  
   limma_results$hit <-
    with(limma_results, ifelse(fdr <= fdr_hit_threshold & abs(logFC) >= log2(fc_hit_threshold), TRUE, FALSE))
  limma_results$hit_annotation <- with(limma_results, 
                                       ifelse(fdr <= fdr_hit_threshold & abs(logFC) >= log2(fc_hit_threshold), 
                                              "hit", 
                                              ifelse(fdr <= fdr_candidate_threshold & abs(logFC) >= log2(fc_candidate_threshold),
                                                     "candidate", "no hit")))
  limma_results$hit_annotation <- factor(limma_results$hit_annotation, ordered = TRUE, levels = c("hit", "candidate", "no hit"))
  return(limma_results)
}

All genetic alterations limma

Conduct limma

limma_results <- NULL

mut_to_limma <- colnames(proteomics_tbl_meta_biomart)[-(1:11)][c(-15)]

for(i in 1:length(mut_to_limma)){
  m <- mut_to_limma[i]
  print(m)
  limma_results <- bind_rows(limma_results, 
                             Calculate_Limma(assay_data = assay(multiomics_MAE[prot_few_nas ,pat_overlap_prot_RNA ,"proteomics"]), mut = m))
}
## [1] "chrom_abber_del11q"
## [1] "chrom_abber_del13q14"
## [1] "chrom_abber_del17p13"
## [1] "chrom_abber_gain8q24"
## [1] "chrom_abber_trisomy12"
## [1] "SNPs_ATM"
## [1] "SNPs_BIRC3"
## [1] "SNPs_EGR2"
## [1] "SNPs_NOTCH1"
## [1] "SNPs_POT1"
## [1] "SNPs_SF3B1"
## [1] "SNPs_TP53"
## [1] "SNPs_XPO1"
## [1] "health_record_bin_IGHV_mutated"
## [1] "health_record_bin_treated"
limma_results$mut %>% unique
##  [1] "chrom_abber_del11q"             "chrom_abber_del13q14"          
##  [3] "chrom_abber_del17p13"           "chrom_abber_gain8q24"          
##  [5] "chrom_abber_trisomy12"          "SNPs_ATM"                      
##  [7] "SNPs_BIRC3"                     "SNPs_EGR2"                     
##  [9] "SNPs_NOTCH1"                    "SNPs_POT1"                     
## [11] "SNPs_SF3B1"                     "SNPs_TP53"                     
## [13] "SNPs_XPO1"                      "health_record_bin_IGHV_mutated"
## [15] "health_record_bin_treated"

Plot limma results

ggplot(data = limma_results) +
  geom_histogram(aes(pvalue.limma, alpha = 0.5), bins = 40) +
  guides(alpha = FALSE) +
  xlab("p-value") +
  facet_wrap( ~ mut +comparison, scale = "free_y") +
  coord_cartesian(xlim = c(0, 1)) +
  pp_sra + ggtitle("p-value histograms")

ggplot(data = limma_results) +
  geom_histogram(aes(pvalue.limma_beforeDEQMS, alpha = 0.5), bins = 40) +
  guides(alpha = FALSE) +
  xlab("p-value") +
  facet_wrap( ~ mut +comparison, scale = "free_y") +
  coord_cartesian(xlim = c(0, 1)) +
  pp_sra + ggtitle("p-values before DEqMS")

Limma Annotation of hits

Annotate

limma_results <- Annotate_Limma_Results(limma_results)

with(limma_results, table(mut, hit_annotation))
##                                 hit_annotation
## mut                               hit candidate no hit
##   chrom_abber_del11q                2         6   7303
##   chrom_abber_del13q14              0         3   7308
##   chrom_abber_del17p13              9        14   7288
##   chrom_abber_gain8q24              1         0   7310
##   chrom_abber_trisomy12           205       332   6774
##   health_record_bin_IGHV_mutated   48       211   7052
##   health_record_bin_treated         3        72   7236
##   SNPs_ATM                          0         0   7311
##   SNPs_BIRC3                        0         4   7307
##   SNPs_EGR2                         4        10   7297
##   SNPs_NOTCH1                       0         0   7311
##   SNPs_POT1                         0         0   7311
##   SNPs_SF3B1                       61        43   7207
##   SNPs_TP53                         4        25   7282
##   SNPs_XPO1                         2         1   7308
message("Number of up and downregulated hits")
## Number of up and downregulated hits
with(limma_results %>% filter(hit==TRUE), table(mut, sign(logFC)))
##                                 
## mut                               -1   1
##   chrom_abber_del11q               2   0
##   chrom_abber_del17p13             8   1
##   chrom_abber_gain8q24             1   0
##   chrom_abber_trisomy12           48 157
##   health_record_bin_IGHV_mutated  21  27
##   health_record_bin_treated        3   0
##   SNPs_EGR2                        0   4
##   SNPs_SF3B1                      53   8
##   SNPs_TP53                        0   4
##   SNPs_XPO1                        2   0
with(limma_results, table(mut, hit)) %>% as_tibble() %>% 
  filter(mut!="health_record_bin_treated") %>%
  mutate(mut=gsub("chrom_abber", "CNVs", mut), mut=gsub("health_record_bin", "IGHV", mut)) %>% 
  separate(mut, into = c("Type", "mut"), sep="_", extra = "merge" ) %>%
  filter(hit==TRUE) %>%
  arrange(desc(n)) %>%
  mutate(mut = as.factor(mut)) %>%
  mutate(mut = factor(mut, levels = .$mut)) %>%
  ggplot(aes(mut, n)) +
  geom_col() +
  ylab("Number of differentially abundant proteins") +
  facet_grid(~Type, scales = "free_x", space="free") +
  pp_sra_noguides_tilted

dif_proteins_P_plot <- with(limma_results %>% mutate(dir=sign(logFC)), table("mut"=paste0(mut,"dir", dir), hit)) %>% as_tibble() %>% 
  separate(mut, into = c("mut", "dir"), sep = "dir") %>%
  filter(mut!="health_record_bin_treated") %>%
  mutate(dir=as.numeric(dir)) %>%
  mutate(n=n*dir) %>%
  mutate(mut=gsub("chrom_abber", "CNVs", mut), mut=gsub("health_record_bin", "IGHV", mut)) %>% 
  separate(mut, into = c("Type", "mut"), sep="_", extra = "merge" ) %>%
  filter(hit==TRUE) %>%
  arrange(desc(n)) %>%
  mutate(mut = as.factor(mut)) %>%
  mutate(mut = factor(mut, levels = unique(.$mut))) %>%
  ggplot(aes(mut, n)) +
  geom_col(aes(fill=as.character(dir))) +
  ylab("Nr. differentially abundant proteins") +
  facet_grid(~Type, scales = "free_x", space="free") +
  pp_sra_noguides_tilted +
  scale_fill_manual(values = c("#0571b0", "#ca0020"), drop=FALSE) +
  geom_hline(yintercept = 0, color="darkgray")

dif_proteins_P_plot

Plot limma hits in volcano plot

All conditions

ggplot(data = limma_results, aes(logFC, -log10(pvalue), colour = hit_annotation)) +
  geom_vline(aes(xintercept = 0)) +
  geom_point() +
  geom_text(aes(label = rowname), 
            data = subset(limma_results, hit_annotation == "hit"), 
            vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
  facet_wrap( ~ mut , ncol = 2) +
  xlab("log2(fold change)") +
  scale_colour_hue() +
  pp_sra

tmp <- limma_results %>% filter(hit==TRUE) %>% 
  mutate(direction=sign(logFC)) %>%
  group_by(mut) %>% arrange(fdr) %>% 
  dplyr::select(rowname, direction)
## Adding missing grouping variables: `mut`
message("Upregulated hits")
## Upregulated hits
sapply(group_nest(tmp)$mut, function(m){
  tmp %>% filter(mut==m, direction==1) %>% .$rowname
})
## $chrom_abber_del11q
## character(0)
## 
## $chrom_abber_del17p13
## [1] "REXO4"
## 
## $chrom_abber_gain8q24
## character(0)
## 
## $chrom_abber_trisomy12
##   [1] "SCYL2"    "EEA1"     "GOLGA3"   "DERA"     "TFCP2"    "EIF4B"   
##   [7] "PTPN11"   "CAMKK2"   "CS"       "TIGAR"    "ANKLE2"   "ESYT1"   
##  [13] "RAB5B"    "ARID2"    "EIF2B1"   "ITPRIP"   "GCN1"     "TBC1D15" 
##  [19] "BRAP"     "DENR"     "USP15"    "PSMD9"    "PTPN6"    "LDHB"    
##  [25] "IRAK4"    "HIP1R"    "TRIM25"   "PYCARD"   "OSBPL8"   "TSFM"    
##  [31] "NAP1L1"   "STRAP"    "RAP1B"    "PPP1R12A" "PYM1"     "PBRM1"   
##  [37] "PA2G4"    "NUDT1"    "PRICKLE1" "PNP"      "C12orf10" "MLX"     
##  [43] "PGAM5"    "SAMHD1"   "GNB4"     "NR2C1"    "SLC2A6"   "RPAP3"   
##  [49] "CUL3"     "SHMT2"    "OGFOD2"   "MON2"     "RILPL2"   "PIKFYVE" 
##  [55] "BTG1"     "CORO1C"   "RAN"      "SART3"    "TRAFD1"   "PUS1"    
##  [61] "TWF1"     "WNK1"     "ARHGDIB"  "CAND1"    "RHOF"     "XPOT"    
##  [67] "USP5"     "MCTS1"    "FRYL"     "HECTD4"   "TXNRD1"   "CCDC91"  
##  [73] "STAT2"    "RAB35"    "ALDH16A1" "AGAP2"    "TAOK3"    "PARP12"  
##  [79] "APPL2"    "IKBIP"    "BRD7"     "SPAG1"    "SCAF11"   "GRPEL1"  
##  [85] "ULK1"     "CKLF"     "PPP2R5B"  "NT5E"     "BCAT1"    "CLIP2"   
##  [91] "ITGAL"    "ACAD10"   "FRMD8"    "SAMD9"    "ATXN2"    "UBE2H"   
##  [97] "ABCD2"    "MADD"     "NLRC4"    "IPO13"    "SCPEP1"   "MLXIP"   
## [103] "SLC25A3"  "ASPSCR1"  "BCL7A"    "ERP29"    "TTC39C"   "AKT2"    
## [109] "ATP6V0A2" "METAP2"   "FES"      "CCDC92"   "NCOR2"    "SETD1B"  
## [115] "RAB21"    "CYFIP1"   "FLNA"     "CNOT8"    "RASSF3"   "LDB1"    
## [121] "TMEM19"   "CALCOCO1" "ANKRD13A" "ITGB2"    "DPYSL2"   "DNAJC7"  
## [127] "ZDHHC17"  "PIP4K2C"  "PRKAB1"   "TMEM231"  "KNTC1"    "EFHD2"   
## [133] "RAB29"    "RHBDF2"   "INF2"     "FIBP"     "PIH1D1"   "C2CD5"   
## [139] "STAT6"    "METTL25"  "RNF135"   "MAP3K8"   "NCKAP5L"  "BIN2"    
## [145] "ITGAX"    "SZRD1"    "TMEM102"  "STK38"    "WASHC3"   "ATF7"    
## [151] "PITPNA"   "UNG"      "HCFC2"    "EEF1A1"   "SH3GL1"   "YARS2"   
## [157] "ATP2A2"  
## 
## $health_record_bin_IGHV_mutated
##  [1] "FTH1"    "FTL"     "USP6NL"  "GRK3"    "TBC1D2B" "ZBTB20"  "ABI3"   
##  [8] "FCRL1"   "SLAMF1"  "FCRL2"   "UBE2E2"  "ADD3"    "DUSP22"  "ADD1"   
## [15] "MTSS1"   "SMC6"    "GNA13"   "KLHL14"  "ELMO2"   "PLEKHF2" "SMC5"   
## [22] "FCHSD2"  "PIP4K2A" "RASA2"   "FLYWCH2" "STK17B"  "POT1"   
## 
## $health_record_bin_treated
## character(0)
## 
## $SNPs_EGR2
## [1] "SLC2A1" "ACTN1"  "EPPK1"  "ZNF446"
## 
## $SNPs_SF3B1
## [1] "SUGP1"    "PSPC1"    "SMARCAD1" "DGCR8"    "ATAD3B"   "EYA3"     "SLC4A1AP"
## [8] "MSH2"    
## 
## $SNPs_TP53
## [1] "LRPAP1" "TP53"   "REXO4"  "FUT11" 
## 
## $SNPs_XPO1
## character(0)
message("Downregulated hits")
## Downregulated hits
sapply(group_nest(tmp)$mut, function(m){
  tmp %>% filter(mut==m, direction==-1) %>% .$rowname
})
## $chrom_abber_del11q
## [1] "CUL5"     "AASDHPPT"
## 
## $chrom_abber_del17p13
## [1] "FXR2"     "MAP2K4"   "PAFAH1B1" "GLOD4"    "ZBTB4"    "ELAC2"    "YWHAE"   
## [8] "RABEP1"  
## 
## $chrom_abber_gain8q24
## [1] "MCPH1"
## 
## $chrom_abber_trisomy12
##  [1] "FRA10AC1" "PELI1"    "TGFBR2"   "TCEAL3"   "BCL7C"    "ABHD14B" 
##  [7] "GDPD1"    "NUMA1"    "SRSF5"    "RANBP3"   "ACYP1"    "ZFP64"   
## [13] "PARP2"    "CRAT"     "TCF4"     "TLK1"     "BCOR"     "WDR89"   
## [19] "PTPRA"    "SELENOM"  "TLE4"     "DNTTIP1"  "CAMK2D"   "MORF4L2" 
## [25] "PCBD1"    "ELMSAN1"  "NOL7"     "GPATCH8"  "UBXN7"    "NCOR1"   
## [31] "PLEKHG1"  "PPP4R3A"  "IRS1"     "SAMSN1"   "TCEAL4"   "ROR1"    
## [37] "C3orf38"  "LANCL1"   "CEP44"    "GPD2"     "ANP32E"   "CTH"     
## [43] "MBD2"     "TTC33"    "DBR1"     "TFAP4"    "SLC25A23" "VAT1"    
## 
## $health_record_bin_IGHV_mutated
##  [1] "CPT1A"    "TMEM165"  "GALNT2"   "UBE2H"    "FAM114A1" "SERPINB6"
##  [7] "UBAP2"    "ATP2C1"   "VAMP5"    "SOGA1"    "ZAP70"    "CCDC186" 
## [13] "ATOX1"    "ATP2A2"   "SLC35B2"  "APOL3"    "PDXDC1"   "PLCH1"   
## [19] "MFN1"     "SEPT10"   "APOBEC3G"
## 
## $health_record_bin_treated
## [1] "MTSS1" "ICAM2" "TNIP1"
## 
## $SNPs_EGR2
## character(0)
## 
## $SNPs_SF3B1
##  [1] "TPP2"     "MAP3K7"   "PPP2R5A"  "WASHC5"   "WASHC4"   "VPS8"    
##  [7] "IQGAP1"   "WASHC1"   "NAA16"    "FKBP15"   "DPH5"     "TAB1"    
## [13] "TTI1"     "CCDC88B"  "SEPSECS"  "TRAPPC6B" "SEPT6"    "TTI2"    
## [19] "WASHC3"   "ABCB7"    "DDX58"    "INPPL1"   "EML3"     "SAFB2"   
## [25] "YWHAB"    "WASHC2A"  "C12orf66" "TTC37"    "SKIV2L"   "DIP2A"   
## [31] "PDS5A"    "NPRL2"    "IDUA"     "RECQL"    "UQCC1"    "TNPO3"   
## [37] "CEP135"   "SZT2"     "PPP6R3"   "AP2A2"    "ARIH1"    "WASHC2C" 
## [43] "JMY"      "C10orf76" "LIG3"     "SUN1"     "SCLT1"    "DEPDC5"  
## [49] "MTMR6"    "KPNA1"    "PNKP"     "MROH1"    "SEPT2"   
## 
## $SNPs_TP53
## character(0)
## 
## $SNPs_XPO1
## [1] "XPO1"  "DHX36"

Trisomy 12 colored by chromosome

limma_results %>% filter(mut=="chrom_abber_trisomy12") %>%
  mutate("Affected_region" = if_else( chromosome_name=="12", "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
  mutate(logFC= if_else(logFC > log2(5), log2(5.2), 
                                    if_else(logFC<(-log2(5)), -log2(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname, color=chromosome_name=="12"), 
            data = subset(limma_results %>% filter(mut=="chrom_abber_trisomy12"), hit_annotation == "hit"), 
            vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  guides(color=FALSE)+
  coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue)))) +
  ggtitle("trisomy12")+
  theme( plot.title = element_text(size = 20)) 

tmp <- limma_results %>% filter(mut=="chrom_abber_trisomy12", chromosome_name=="12") %>%
  .$hit_annotation %>% table %>% as_tibble() 
colnames(tmp) <- c("Protein", "n")
tmp
tmp %>%
  mutate(Protein=as.factor(Protein)) %>% 
  mutate(Protein=factor(Protein, levels = c("hit", "candidate", "no hit"))) %>%
  ggplot(aes(x="", fill=Protein,y= n)) +
  geom_col() +
  ggtitle("Hit annotation proteins on chromosome 12 in trisomy 12") + 
  coord_polar("y", start=0) +
  pp_sra +
  theme(axis.title = element_blank(), axis.ticks = element_blank())

tmp <- (limma_results %>% 
              filter(mut=="chrom_abber_trisomy12", chromosome_name=="12") %>%
  .$logFC>0)  %>% table %>% as_tibble() 
colnames(tmp) <- c("fc", "n")
tmp <- tmp %>% arrange(desc(`fc`)) 
tmp <- tmp %>% mutate(label=paste( round(n/sum(tmp$n)*100, 1), "%" ), cums =cumsum(tmp$n) ) 
tmp %>%
  mutate("fc"=as.factor(`fc`)) %>% 
  ggplot(aes(x="", fill=`fc`,y= n)) +
  geom_col() +
  ggtitle("Proteins on chromosome 12 with positive logFC in trisomy 12") + 
  coord_polar("y", start=0) +
  pp_sra +
  theme(axis.title = element_blank(), axis.ticks = element_blank())+
  scale_fill_manual(values = c( "grey", "#0571b0"), drop=FALSE)+
  annotate(geom = "text", y = tmp$cums-(tmp$n/2), x = 1, label = tmp$label)

tmp <- (limma_results %>%  
  filter(mut=="chrom_abber_trisomy12", hit==TRUE) %>%
  .$chromosome_name == "12") %>%
  table %>% as_tibble()
colnames(tmp) <- c("chr12", "n")
tmp <- tmp %>% arrange(desc(`chr12`)) 
tmp <- tmp %>% mutate(label=paste( round(n/sum(tmp$n)*100, 1), "%" ), cums =cumsum(tmp$n) ) 
piechart_hits_prot_chr_tris12 <-  tmp %>%
  mutate("chr12"=as.factor(`chr12`)) %>% 
  ggplot(aes(x="", fill=`chr12`,y= n)) +
  geom_col() +
  ggtitle("Chromosome location of hits in trisomy 12") + 
  coord_polar("y", start=0) +
  pp_sra +
  theme(axis.title = element_blank(), axis.ticks = element_blank())+
  scale_fill_manual(values = c( "grey", "#0571b0"), drop=FALSE)
piechart_hits_prot_chr_tris12 +
  annotate(geom = "text", y = tmp$cums-(tmp$n/2), x = 1, label = tmp$label)

del13q14 colored by chromosome

limma_results %>% filter(mut=="chrom_abber_del13q14") %>%
  mutate("Affected_region" = if_else( (chromosome_name=="13" & start_position>47000000 & start_position<51000000 ), "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
  mutate(logFC= if_else(logFC > log2(5), log2(5.2), 
                                    if_else(logFC<(-log2(5)), -log(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname, color=(chromosome_name=="13" & start_position>47000000 & start_position<51000000)), 
            data = subset(limma_results %>% filter(mut=="chrom_abber_del13q14"), hit_annotation == "hit"), 
            vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  guides(color=FALSE)+
  coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue))))+
  ggtitle("del13q14")+
  theme( plot.title = element_text(size = 20)) 

del11q colored by chromosome

limma_results %>% filter(mut=="chrom_abber_del11q") %>%
  mutate("Affected_region" = if_else( (chromosome_name=="11" & start_position>97400000 & start_position<110600000), "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
  mutate(logFC= if_else(logFC > log2(5), log2(5.2), 
                                    if_else(logFC<(-log2(5)), -log2(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname, color=(chromosome_name=="11" & start_position>97400000 & start_position<110600000)), 
            data = subset(limma_results %>% filter(mut=="chrom_abber_del11q"), hit_annotation == "hit"), 
            vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  guides(color=FALSE)+
  coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue)))) +
  ggtitle("del11q")+
  theme(plot.title = element_text(size = 20)) 

del17p13 colored by chromosome

limma_results %>% filter(mut=="chrom_abber_del17p13") %>%
  mutate("Affected_region" = if_else( (chromosome_name=="17" & mean_position<10800000), "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
  mutate(logFC= if_else(logFC > log2(5), log2(5.2), 
                                    if_else(logFC<(-log2(5)), -log2(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname, color=(chromosome_name=="17" & mean_position<10800000)), 
            data = subset(limma_results %>% filter(mut=="chrom_abber_del17p13"), hit_annotation == "hit"), 
            vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  guides(color=FALSE)+
  coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue)))) +
  ggtitle("del17p13")+
  theme(plot.title = element_text(size = 20)) 

XPO1

XPO_volcano_plot <- limma_results %>% filter(mut=="SNPs_XPO1") %>%
  mutate("Affected_region" = if_else( rowname=="XPO1", "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname), 
            data = subset(limma_results %>% filter(mut=="SNPs_XPO1", rowname!="XPO1"), hit_annotation == "hit"), 
            vjust = 0, size = 3, nudge_y = (-1), nudge_x = -0.3,#  nudge_y = -0.3, 
            check_overlap = FALSE, color="#0571b0") +
   geom_text(aes(label = rowname, color=rowname=="XPO1"), 
            data = subset(limma_results %>% filter(mut=="SNPs_XPO1", rowname=="XPO1")), 
            vjust = 0, size = 3,  nudge_x = 0.2, nudge_y = (-0.5), #nudge_y = (-0.1), 
            check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  coord_cartesian(xlim=c(-log2(5), log2(5)) ) 

XPO_volcano_plot + ggtitle("XPO1")+
  theme( plot.title = element_text(size = 20)) +
  guides(color=FALSE)

TP53

TP53_volcano_plot <- limma_results %>% filter(mut=="SNPs_TP53") %>%
  mutate("Affected_region" = if_else( rowname=="TP53", "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname), 
            data = subset(limma_results %>% filter(mut=="SNPs_TP53", rowname!="TP53"), hit_annotation == "hit"), 
            vjust = 0, size = 3,  nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE, color="#0571b0") +
   geom_text(aes(label = rowname, color=rowname=="TP53"), 
            data = subset(limma_results %>% filter(mut=="SNPs_TP53", rowname=="TP53")), 
            vjust = 0, size = 3, nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  coord_cartesian(xlim=c(-log2(2.2), log2(2.2)) ) 

TP53_volcano_plot + ggtitle("TP53")+
  theme( plot.title = element_text(size = 20)) +
  guides(color=FALSE)

ATM

ATM_volcano_plot <- limma_results %>% filter(mut=="SNPs_ATM") %>%
  mutate("Affected_region" = if_else( rowname=="ATM", "TRUE", 
                                     if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname), 
            data = subset(limma_results %>% filter(mut=="SNPs_ATM", rowname!="ATM"), hit_annotation == "hit"), 
            vjust = 0, size = 3,  nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE, color="#0571b0") +
   geom_text(aes(label = rowname, color=rowname=="ATM"), 
            data = subset(limma_results %>% filter(mut=="SNPs_ATM", rowname=="ATM")), 
            vjust = 0, size = 3, nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
  pp_sra+
  coord_cartesian(xlim=c(-log2(2.2), log2(2.2)) ) 

ATM_volcano_plot + ggtitle("ATM")+
  theme( plot.title = element_text(size = 20)) +
  guides(color=FALSE)

Limma with subsets of patients

Consensus clusters groups 5 vs. all others

selpats <- colData(multiomics_MAE) %>% as_tibble() %>% filter(!is.na(PG)) %>% .$patient_ID
proteomics_tbl_meta_biomart_tmp <- proteomics_tbl_meta_biomart
proteomics_tbl_meta_biomart <- left_join(proteomics_tbl_meta_biomart, enframe(CCP_group5==5, name = "primary", value="PG"), by="primary" )
limma_results_CC5_5_all <-  Calculate_Limma(assay_data = assay(multiomics_MAE[prot_few_nas , selpats ,"proteomics"]), mut = "PG")
## ExperimentList contains data.frame or DataFrame,
##   potential for errors with mixed data types
## ExperimentList contains data.frame or DataFrame,
##   potential for errors with mixed data types
proteomics_tbl_meta_biomart <- proteomics_tbl_meta_biomart_tmp 

ggplot(data = limma_results_CC5_5_all) +
  geom_histogram(aes(pvalue.limma, alpha = 0.5), bins = 40) +
  guides(alpha = FALSE) +
  xlab("p-value") +
  facet_wrap( ~ mut , scale = "free_y") +
  coord_cartesian(xlim = c(0, 1)) +
  pp_sra +
  ggtitle("PG groups")

limma_results_CC5_5_all <- Annotate_Limma_Results(limma_results_CC5_5_all)

limma_results_CC5_5_all %>% 
  mutate("Affected_region" =  if_else(hit_annotation == "hit",  "hit", "not altered")) %>% 
  mutate("Affected_region"= factor(Affected_region, levels=c("hit", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue)))  +
  geom_vline(aes(xintercept = 0)) +
  geom_point(alpha=0.8, aes(colour = Affected_region)) +
  geom_text(aes(label = rowname), color="#0571b0", 
            data = subset(limma_results_CC5_5_all %>% filter( hit_annotation == "hit"), 
            vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE)) +
  facet_wrap( ~ mut, ncol = 1) +
  xlab("log2(fold change)") +
  scale_color_manual(values = c("#0571b0",  "grey"), drop=FALSE) +
  pp_sra+
  guides(color=guide_legend(title = "PG groups"))+
  ggtitle("Consensus cluster plus group 5 against all others") 

limma_results_CC5_5_all %>% filter(hit==TRUE) %>% .$rowname
##   [1] "WDR1"      "ACTR2"     "ACTR3"     "FNTA"      "CAP1"      "ATG5"     
##   [7] "SH3BGRL"   "GAPVD1"    "NRBP1"     "CDC5L"     "PITPNB"    "BTK"      
##  [13] "ARPC3"     "TERF2"     "ARPC1B"    "SARS"      "LUC7L3"    "ATG7"     
##  [19] "VAV2"      "CFL1"      "ARPC5"     "IKBKB"     "NUP155"    "PGK1"     
##  [25] "TBC1D7"    "CWC25"     "DHX9"      "PPP6C"     "PCF11"     "PANK4"    
##  [31] "PIK3CD"    "NUP98"     "TERF2IP"   "AFP"       "NUP54"     "SRSF4"    
##  [37] "ARPC2"     "RALGAPB"   "PDP1"      "NUP93"     "PUF60"     "OTULIN"   
##  [43] "GLE1"      "CSK"       "GRK2"      "NUP133"    "RANBP1"    "ARHGDIA"  
##  [49] "USP9X"     "ARHGAP17"  "ANKRD44"   "SNX5"      "SF3B1"     "RABGGTB"  
##  [55] "UBR1"      "C16orf58"  "HNRNPC"    "ZNRF2"     "HTT"       "ARFGEF1"  
##  [61] "SF3B6"     "PNN"       "CRNKL1"    "SFPQ"      "GSR"       "CORO1A"   
##  [67] "NUP58"     "PIK3R1"    "RBM12B"    "METAP2"    "TWF2"      "SH3GL1"   
##  [73] "OSTF1"     "NUP205"    "MOB2"      "YWHAB"     "TLN1"      "PPP6R1"   
##  [79] "SLTM"      "OTUB1"     "EMC3"      "PAXBP1"    "HNRNPR"    "SMPD4"    
##  [85] "STAT6"     "MTX1"      "PGGT1B"    "PRPF19"    "DDX46"     "WIPF1"    
##  [91] "TRAPPC10"  "USP47"     "PRKAG1"    "CCDC47"    "NUP210"    "MAPRE1"   
##  [97] "EEF2"      "ASXL2"     "RPN2"      "BUD13"     "RAC2"      "UBE2O"    
## [103] "CRKL"      "CAPZB"     "GDI2"      "MTMR9"     "RPN1"      "SNRNP48"  
## [109] "PAF1"      "KDELC1"    "NXF1"      "THRAP3"    "HECTD1"    "SON"      
## [115] "CIAO1"     "EEF1A1"    "TBCA"      "RASSF2"    "AKT2"      "ATG3"     
## [121] "SAP18"     "MAP2K2"    "SERPINC1"  "UBE2V2"    "DENND1C"   "YEATS2"   
## [127] "NONO"      "EFHD2"     "FBXW2"     "NDUFAF2"   "PAPOLG"    "YLPM1"    
## [133] "SUPT6H"    "SSR4"      "MARK2"     "UBE2N"     "TNFAIP8"   "BIN2"     
## [139] "HAX1"      "PTGES3"    "SRSF3"     "DIAPH1"    "SF3B2"     "IMMT"     
## [145] "IST1"      "UBAC1"     "ZFR"       "DCTN4"     "MTMR6"     "YKT6"     
## [151] "ZNF512"    "CPSF3"     "KHDRBS1"   "PTPA"      "GTPBP1"    "CDC73"    
## [157] "PDPR"      "TF"        "TSC1"      "CCDC25"    "WASF2"     "BRK1"     
## [163] "HMGXB4"    "GRB2"      "RSRC1"     "CLEC3B"    "MYD88"     "PLEKHM2"  
## [169] "INPP5D"    "NUP107"    "HNRNPM"    "ZBTB48"    "C16orf70"  "GFM1"     
## [175] "TPD52L2"   "RBBP6"     "PTPN11"    "SNRNP70"   "FBLN1"     "PDCD6IP"  
## [181] "CAPZA1"    "YTHDC1"    "NFKBIE"    "DARS2"     "CYC1"      "DDHD2"    
## [187] "COMMD10"   "PFN1"      "NUP62"     "MYL6"      "PSMB2"     "NCKAP1L"  
## [193] "EMC1"      "STAT5A"    "HM13"      "RNPS1"     "UVRAG"     "PACS1"    
## [199] "ILF3"      "XPNPEP1"   "PIK3C3"    "INPPL1"    "CARS"      "PNPLA8"   
## [205] "NFKB1"     "DDX10"     "DBNL"      "NAP1L4"    "CERS2"     "TOMM70"   
## [211] "SF3A3"     "VPS4B"     "SF3A1"     "RBM17"     "BAP1"      "GNPDA2"   
## [217] "TRAPPC1"   "EXOC4"     "ARSA"      "CAMLG"     "GSPT1"     "POLE4"    
## [223] "GLS"       "MBIP"      "NUP88"     "NUPL2"     "LRBA"      "DCTN6"    
## [229] "APMAP"     "CHMP4B"    "UTP3"      "MAP2K1"    "POMGNT1"   "GMFB"     
## [235] "SPCS2"     "LTA4H"     "RDH14"     "PRKAA1"    "RALGAPA2"  "GAK"      
## [241] "CENPB"     "KYAT1"     "RBM26"     "CYFIP2"    "SLC25A12"  "UBA5"     
## [247] "NRAS"      "USP24"     "DDX50"     "PPP3R1"    "IFI16"     "ADAR"     
## [253] "REL"       "ATP6V1C1"  "NUP188"    "PCCB"      "MSL1"      "SRRM2"    
## [259] "RMC1"      "RABL6"     "RIF1"      "ELP2"      "RNF6"      "EMC2"     
## [265] "WDR37"     "SUZ12"     "TMEM214"   "RIC1"      "SRSF11"    "HIBADH"   
## [271] "TRAPPC9"   "PREP"      "PDPK1"     "RANBP2"    "PPP3CB"    "VAPA"     
## [277] "CDK14"     "DOCK11"    "TSC2"      "PRPSAP1"   "OFD1"      "DKC1"     
## [283] "WAS"       "PIBF1"     "OCIAD1"    "EXOC2"     "RABGGTA"   "ALDH5A1"  
## [289] "CTR9"      "TPR"       "HNRNPH2"   "HSH2D"     "METTL21A"  "FAM160B1" 
## [295] "GOLGA5"    "RNF126"    "DOCK2"     "ESF1"      "PLA2G15"   "DCTN2"    
## [301] "RBM39"     "SSR3"      "ACAP2"     "RBM15"     "ATP6V1B2"  "IARS2"    
## [307] "ITIH2"     "CCT7"      "BRAP"      "RPS12"     "DNMT1"     "TCP1"     
## [313] "PLCG2"     "MOSPD2"    "AKAP13"    "ATP6V1H"   "DAPP1"     "TFIP11"   
## [319] "TMX1"      "ZFC3H1"    "CLPB"      "SNRPA1"    "CISD2"     "CNN2"     
## [325] "MALT1"     "MAPK1"     "HABP2"     "ENOPH1"    "DOK1"      "PCCA"     
## [331] "MFN2"      "PTP4A2"    "SLU7"      "CCAR1"     "ANKRD13A"  "PSMA6"    
## [337] "EGLN1"     "FLII"      "COMMD9"    "ELMOD2"    "SAFB"      "PKN1"     
## [343] "FKBP1A"    "CYB5B"     "PRPS1"     "CAND1"     "NOL9"      "PGK2"     
## [349] "TRMT61B"   "NTMT1"     "GGCT"      "VAV1"      "SNRPB2"    "ATP7A"    
## [355] "HCLS1"     "NELFCD"    "HACD3"     "PDLIM2"    "ARAF"      "PTPN6"    
## [361] "PDCD7"     "CLEC16A"   "MLYCD"     "PDK3"      "PES1"      "GPN1"     
## [367] "ARMT1"     "SNW1"      "CBX3"      "SRSF5"     "AKT1"      "BECN1"    
## [373] "SRSF10"    "MAPK14"    "CENPC"     "CPSF2"     "PBDC1"     "TECPR1"   
## [379] "APPL1"     "A2M"       "SLC25A36"  "PLAA"      "AIFM1"     "ZMPSTE24" 
## [385] "ROCK1"     "MAPK15"    "GMFG"      "TMEM208"   "BAZ1B"     "EMC7"     
## [391] "ETF1"      "PTEN"      "TRAM1"     "PDHA1"     "HOOK3"     "RABGAP1"  
## [397] "FUNDC2"    "VAMP4"     "MCM3AP"    "SMARCA5"   "TRADD"     "MSN"      
## [403] "TST"       "IKBKG"     "RMDN1"     "UNC45A"    "RFX1"      "APAF1"    
## [409] "ATG16L1"   "IK"        "ATP6V1A"   "WDR44"     "CLASRP"    "NUP35"    
## [415] "RPL12"     "CSTF3"     "TRIP11"    "ABL1"      "CANX"      "NIPSNAP2" 
## [421] "STIM2"     "UBR2"      "PAWR"      "ISYNA1"    "PLEKHF2"   "PLS1"     
## [427] "SEC22B"    "KMT2B"     "TMCC1"     "RABAC1"    "DDRGK1"    "KIAA0100" 
## [433] "CAPN1"     "ABI1"      "TTC9"      "GNL1"      "TFAM"      "STT3B"    
## [439] "GLTP"      "VPS33A"    "RASGRP2"   "HMGN3"     "MTMR12"    "IRAK4"    
## [445] "ZC3H18"    "PDHB"      "PAAF1"     "DCTN5"     "GAR1"      "SPTY2D1"  
## [451] "DHX16"     "VPS37A"    "GOT1"      "WTAP"      "LIMD1"     "WDR5"     
## [457] "AP3M1"     "ELP4"      "CLMN"      "SPCS3"     "RSBN1L"    "RPL14"    
## [463] "RHOT1"     "SLC25A32"  "ITCH"      "DDOST"     "IKZF1"     "DNAJC11"  
## [469] "CHCHD6"    "SFXN3"     "TNIP2"     "CLIC1"     "EML3"      "ARHGDIB"  
## [475] "ERGIC1"    "EXOC3"     "CEP44"     "AUP1"      "EXOC1"     "RNF123"   
## [481] "MBLAC2"    "CAPNS1"    "RPE"       "ALDH18A1"  "HTATIP2"   "RAD23A"   
## [487] "VDAC1"     "VCPIP1"    "BPTF"      "DENND4C"   "GPATCH4"   "SF3B4"    
## [493] "CCDC90B"   "UBAP1"     "IWS1"      "RBM10"     "CTNNBL1"   "TARS"     
## [499] "NDC1"      "PIK3R4"    "DRG2"      "DIDO1"     "CCT8"      "SLC38A2"  
## [505] "U2SURP"    "BCL10"     "RIC8A"     "SF1"       "MTPN"      "BLOC1S3"  
## [511] "MAPKAPK5"  "GATAD1"    "MED1"      "SRSF7"     "PTK2B"     "MUL1"     
## [517] "EIF5"      "GSTCD"     "MTMR14"    "PLRG1"     "NECAP2"    "STAT3"    
## [523] "XRN2"      "BLOC1S5"   "PIGX"      "SSRP1"     "CHCHD3"    "EMC4"     
## [529] "CHD4"      "LSG1"      "RACK1"     "FBXO21"    "SYMPK"     "PPIA"     
## [535] "CUL3"      "ZNF638"    "FAM76A"    "SH2D3C"    "COTL1"     "TBCB"     
## [541] "HADHB"     "MORC2"     "BRAF"      "CHMP2B"    "ATP6AP2"   "SMC3"     
## [547] "XAB2"      "CLPTM1L"   "MPHOSPH10" "NELFA"     "KYNU"      "SRSF1"    
## [553] "EXOC7"     "LMAN2"     "MTIF3"     "NFATC2"    "PSMB9"     "RAD21"    
## [559] "AKT3"      "SEC61B"    "MCU"       "MTDH"      "SRRM1"     "YWHAG"    
## [565] "FGD3"      "ELMO2"     "BLOC1S2"   "C2CD2L"    "TRAPPC6B"  "OGFR"     
## [571] "PDZD8"     "MAD1L1"    "ILF2"      "RALGPS2"   "ELP3"      "TXNRD1"   
## [577] "EPM2A"     "UBE2Z"     "MYL12A"    "OLA1"      "FUNDC1"    "RETSAT"   
## [583] "PREB"      "CDC42EP3"  "STXBP2"    "RPRD2"     "BLOC1S4"   "DCTN3"    
## [589] "LRRFIP1"   "CDK7"      "SNRPD2"    "LCP1"      "CLPX"      "ATG4B"    
## [595] "CDC42"     "AKAP1"     "GEMIN5"    "GFM2"      "KIAA0753"  "PFKL"     
## [601] "PI4KB"     "CAB39"     "CCT3"      "NKIRAS2"   "DCTN1"     "GON4L"    
## [607] "UBE2M"     "RABGEF1"   "TRMT1L"    "DUT"       "MRPL15"    "NAB1"     
## [613] "UBA6"      "RGL2"      "GPATCH8"   "ANKMY2"    "DTNBP1"    "PPP1R12C" 
## [619] "CAPN7"     "RPL30"     "VMP1"      "COPS7B"    "GATAD2B"   "TRAPPC3"  
## [625] "ZC3H4"     "FERMT3"    "SIN3B"     "PPP3CC"    "RPSA"      "KPNA2"    
## [631] "ZCCHC6"    "MYO1E"     "GOLIM4"    "ATAD2"     "VPS35L"    "HADH"     
## [637] "KIF3B"     "SPOP"      "BIN1"      "PRKAB1"    "PRKD2"     "RPL22"    
## [643] "TMED10"    "PSMD8"     "MRGBP"     "PRKCZ"     "RRAS2"     "ATP6V1G1" 
## [649] "MSANTD4"   "AHCTF1"    "YTHDF3"    "ARPC5L"    "UBR4"      "PTPN7"    
## [655] "CEP162"    "SNX2"      "CBFB"      "NARS"      "STT3A"     "NSDHL"    
## [661] "CDKN2AIP"  "KIF5B"     "ACOX1"     "SLC22A18"  "ARIH1"     "NIPSNAP1" 
## [667] "MCUB"      "NAP1L1"    "RP2"       "KIFAP3"    "SREK1"     "TPP1"     
## [673] "UFL1"      "PPOX"      "PPP1R16B"  "STAP1"     "ELP1"      "CNOT9"    
## [679] "NDUFAF6"   "TCERG1"    "VIRMA"     "TARS2"     "MDH2"      "GRK3"     
## [685] "NUP160"    "SERPINF2"  "USP25"     "CLIC4"     "RBMX"      "RRP9"     
## [691] "RRP1"      "PRDX3"     "IMP3"      "AIMP2"     "BTF3L4"    "PHPT1"    
## [697] "TBC1D9B"   "CNDP2"     "SSBP1"     "WASHC3"    "PRAG1"     "FMC1"     
## [703] "TTI1"      "PDHX"      "TOM1"      "NSMAF"     "NAPEPLD"   "GDI1"     
## [709] "MAPK9"     "ZNF326"    "MICU1"     "RIPOR2"    "NCOA5"     "NLN"      
## [715] "SNX3"      "LYAR"      "ZBTB9"     "NUP50"     "FNBP4"     "C3orf33"  
## [721] "GPN2"      "DOPEY2"    "PSMA5"     "ADRM1"     "FAF2"      "SH3BP2"   
## [727] "ZNHIT1"    "MIEN1"     "CHRAC1"    "SMARCE1"   "CCDC43"    "PIGK"     
## [733] "ERCC4"     "MOB3A"     "KNOP1"     "RPL38"     "EXOC8"     "PXK"      
## [739] "XXYLT1"    "SLC4A1AP"  "HMGCL"     "MTX2"      "LIMD2"     "TMED2"    
## [745] "MAN1B1"    "MYO9B"     "HNRNPH3"   "MAP3K8"    "USP3"      "RIPK1"    
## [751] "IGBP1"     "GNB2"      "ATP6V1E1"  "CIRBP"

Gene set enrichment analysis consensus cluster plus groups 5 vs. all

####### Input of all Values
longFormat((multiomics_MAE[prot_few_nas , selpats ,"proteomics"])) %>% 
  as_tibble() %>% 
  dplyr::select("NAME"=rowname, "DESCRIPTION"=assay , primary, value) %>%
  unique() %>%
  spread(primary, value) #%>%
  #write.table(file = "/Users/sophierabe/Desktop/PhD/Labor/Proteomics/CLL/GSEA_CLL_Proteomics/190927_GSEA_data_long_gradient_CCP_5vsall.txt", quote = FALSE, row.names = FALSE, sep = "\t", na="")

tmp_pats <- longFormat((multiomics_MAE[prot_few_nas , selpats ,"proteomics"])) %>% 
  as_tibble() %>% 
  dplyr::select("NAME"=rowname, primary, value) %>%
  unique() %>%
  spread(primary, value)  %>%
  dplyr::select(-1) %>%
  colnames() %>%
  enframe() 
tmp_pats <- tmp_pats %>%
  left_join(., 
            enframe(multiomics_MAE$PG, value="PG", name="primary" ) %>% mutate("PG"=PG==5) , 
            by=c("value"="primary") )
tmp_pats %>%
  dplyr::select(PG) #%>% t() %>%
  #write.table(file = "/Users/sophierabe/Desktop/PhD/Labor/Proteomics/CLL/GSEA_CLL_Proteomics/190927_GSEA_phenotypelabels_long_gradient__CCP_5vsall.cls", quote = FALSE, row.names = FALSE, col.names = FALSE, sep = " ", na="")
dim(tmp_pats)
sum(tmp_pats$PG==TRUE)

Save data needed for further analysis

save(limma_results,
     file = "/Volumes/sd17b003/Sophie/Analysis/CLL_Proteomics/CLL_Proteomics_final/Proteomics_Git/Robjects/CLL_Proteomics_LimmaProteomics.RData")

Session Info

sessionInfo()
## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
## 
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] DEqMS_1.0.1                 cowplot_1.0.0              
##  [3] survminer_0.4.6             ConsensusClusterPlus_1.46.0
##  [5] RColorBrewer_1.1-2          gplots_3.0.1.1             
##  [7] PMA_1.1                     MultiAssayExperiment_1.8.3 
##  [9] biomaRt_2.38.0              biomartr_0.9.0             
## [11] Rtsne_0.15                  IHW_1.10.1                 
## [13] apeglm_1.4.2                limma_3.38.3               
## [15] fdrtool_1.2.15              vsn_3.50.0                 
## [17] forcats_0.4.0               stringr_1.4.0              
## [19] dplyr_0.8.3                 purrr_0.3.3                
## [21] readr_1.3.1                 tidyr_1.0.0                
## [23] tibble_2.1.3                tidyverse_1.3.0            
## [25] TOSTER_0.3.4                DESeq2_1.22.2              
## [27] SummarizedExperiment_1.12.0 DelayedArray_0.8.0         
## [29] BiocParallel_1.16.6         matrixStats_0.55.0         
## [31] Biobase_2.42.0              GenomicRanges_1.34.0       
## [33] GenomeInfoDb_1.18.2         IRanges_2.16.0             
## [35] S4Vectors_0.20.1            BiocGenerics_0.28.0        
## [37] readxl_1.3.1                ggrepel_0.8.1              
## [39] ggpmisc_0.3.3               ggfortify_0.4.8            
## [41] ggbeeswarm_0.6.0            effsize_0.7.6              
## [43] ggpubr_0.2.4                magrittr_1.5               
## [45] lemon_0.4.3                 Hmisc_4.3-0                
## [47] ggplot2_3.2.1               Formula_1.2-3              
## [49] survival_3.1-8              lattice_0.20-38            
## [51] Matrix_1.2-18               progress_1.2.2             
## [53] reshape2_1.4.3              pheatmap_1.0.12            
## [55] openxlsx_4.1.4              gtools_3.8.1               
## [57] plyr_1.8.4                 
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_0.2.5       RSQLite_2.1.4          AnnotationDbi_1.44.0  
##   [4] htmlwidgets_1.5.1      munsell_0.5.0          codetools_0.2-16      
##   [7] preprocessCore_1.44.0  withr_2.1.2            colorspace_1.4-1      
##  [10] knitr_1.26             rstudioapi_0.10        ggsignif_0.6.0        
##  [13] labeling_0.3           slam_0.1-46            bbmle_1.0.20          
##  [16] GenomeInfoDbData_1.2.0 lpsymphony_1.10.0      KMsurv_0.1-5          
##  [19] farver_2.0.1           bit64_0.9-7            coda_0.19-3           
##  [22] vctrs_0.2.0            generics_0.0.2         xfun_0.11             
##  [25] R6_2.4.1               locfit_1.5-9.1         bitops_1.0-6          
##  [28] assertthat_0.2.1       scales_1.1.0           nnet_7.3-12           
##  [31] beeswarm_0.2.3         gtable_0.3.0           affy_1.60.0           
##  [34] rlang_0.4.2            zeallot_0.1.0          genefilter_1.64.0     
##  [37] splines_3.5.2          lazyeval_0.2.2         acepack_1.4.1         
##  [40] impute_1.56.0          broom_0.5.2            checkmate_1.9.4       
##  [43] BiocManager_1.30.10    yaml_2.2.0             modelr_0.1.5          
##  [46] backports_1.1.5        tools_3.5.2            ellipsis_0.3.0        
##  [49] affyio_1.52.0          Rcpp_1.0.3             base64enc_0.1-3       
##  [52] zlibbioc_1.28.0        RCurl_1.95-4.12        prettyunits_1.0.2     
##  [55] rpart_4.1-15           zoo_1.8-6              haven_2.2.0           
##  [58] cluster_2.1.0          fs_1.3.1               data.table_1.12.8     
##  [61] reprex_0.3.0           hms_0.5.2              evaluate_0.14         
##  [64] xtable_1.8-4           XML_3.98-1.20          emdbook_1.3.11        
##  [67] gridExtra_2.3          compiler_3.5.2         KernSmooth_2.23-16    
##  [70] crayon_1.3.4           htmltools_0.4.0        geneplotter_1.60.0    
##  [73] lubridate_1.7.4        DBI_1.0.0              dbplyr_1.4.2          
##  [76] MASS_7.3-51.4          cli_2.0.0              gdata_2.18.0          
##  [79] pkgconfig_2.0.3        km.ci_0.5-2            numDeriv_2016.8-1.1   
##  [82] foreign_0.8-72         xml2_1.2.2             annotate_1.60.1       
##  [85] vipor_0.4.5            XVector_0.22.0         rvest_0.3.5           
##  [88] digest_0.6.23          Biostrings_2.50.2      rmarkdown_1.18        
##  [91] cellranger_1.1.0       survMisc_0.5.5         htmlTable_1.13.3      
##  [94] curl_4.3               lifecycle_0.1.0        nlme_3.1-142          
##  [97] jsonlite_1.6           fansi_0.4.0            pillar_1.4.2          
## [100] httr_1.4.1             glue_1.3.1             zip_2.0.4             
## [103] bit_1.1-14             stringi_1.4.3          blob_1.2.0            
## [106] latticeExtra_0.6-28    caTools_1.17.1.3       memoise_1.1.0